clc; clear; close all;

% 参数定义
alpha = 2.4;
A = [-alpha, 0, 0;
     0, 0, 0;
     0, 0, 0];
B = [1, -4, -3.5;
     0, 1, 2;
    -1, -4, 1.5];
k = 0.2;

% 固定phi0
phi0 = -0.58;

% 时间参数
dt = 0.01;
T_end = 4000;
N = floor(T_end / dt);

% 微分方程句柄
f = @(t,x) mCNN(t, x, A, B, k);

% 初始条件向量
x0 = [1e-6; 0; 0; 0; 0; 0; 0; 0; phi0];

% 预分配存储矩阵
X = zeros(length(x0), N+1);
X(:,1) = x0;

% RK4积分
for i = 1:N
    X(:, i+1) = RK4(f, X(:, i), dt);
end

% 时间向量
T = linspace(0, T_end, N+1);

% 选取时间区间 3900~4000s
sel = T >= 3900;
X_sel = X(:, sel);
T_sel = T(sel);

% 绘图
figure('Position',[100 100 1200 400]);

title(sprintf('\\phi_0 = %.2f', phi0));

% a1 相空间 x11-x13 和 x21-x23
subplot(1,3,1);
plot(X_sel(1,:), X_sel(3,:), 'r'); hold on;
plot(X_sel(4,:), X_sel(6,:), 'b--');
xlabel('x_{11}, x_{21}');
ylabel('x_{13}, x_{23}');
legend('x_{11}-x_{13}', 'x_{21}-x_{23}');
grid on;

% a2 时间序列 x13 和 x23
subplot(1,3,2);
plot(T_sel, X_sel(3,:), 'r'); hold on;
plot(T_sel, X_sel(6,:), 'b--');
xlabel('t (s)');
ylabel('x_{13}, x_{23}');
legend('x_{13}', 'x_{23}');
grid on;

% a3 相空间 x13-x23
subplot(1,3,3);
plot(X_sel(3,:), X_sel(6,:), 'b');
xlabel('x_{13}');
ylabel('x_{23}');
grid on;